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1 Introduction 

The purpose of measuring ECGs at the body surface 
or MCGs outside the body is to obtain information 
of the functioning of the heart. The electric sources 
within the heart, that generate both the ECG and the 
MCG, cannot be determined uniquely from the mea¬ 
sured signals. This is a result of the fact that there 
is an infinite number of source distributions (equiva¬ 
lent sources) that reproduce the same external fields 
as the true source. As a consequence, a source model 
is needed in order to constrain the source description 
in such a way that there is a unique solution to the 
inverse problem. 

In order to be clinically useful, the source model has 
be physiological realistic and yield source parame¬ 
ters that allow a physiological interpretation. In sit¬ 
uations where the electric activity is confined to a 
small region (such as early activation in the heart) 
a simple source model like a moving current dipole 
is adequate: the source parameters (i.e. dipole po¬ 
sition, strength and direction) have a clear physio¬ 
logical meaning. However, throughout the cardiac 
cycle activity spreads in a complex pattern over the 
myocardium, and although a single moving dipole 
can explain the measured fields in this situation quite 
well, the physiological interpretation of the source pa¬ 
rameters is obscure. 

For the adequate description of the development of 
the electric activity during the cardiac cycle dis¬ 
tributed source models are needed, such as the epicar- 
dial potential distribution [1] or the Uniform Double 
Layer (UDL) model [2], which was used in this study. 
In the UDL model the cardiac source during ventricu¬ 
lar depolarization (the QRS complex) is described by 
the local activation times at the ventricular surface. 
Distributed models have many degrees of freedom 
(i.e. a large number of source parameters). In a noise- 
free world, there is a unique relation between these 
source parameters and the external fields. However, 
many different sets of source parameters generate ex¬ 


ternal fields that differ less than realistic noise levels. 
As a consequence the inverse problem involved is ill- 
posed, and further constraints on the source param¬ 
eters are needed, effectively reducing the number of 
degrees of freedom. Commonly these constraints are 
implemented as regularization: the extend to which 
the source parameters do not comply to some desir¬ 
able property is expressed in a regularization func¬ 
tion, and rather than just minimizing the difference 
between the measured and model fields, the weighted 
sum of this difference and the regularization function 
is minimized. 

The extend to which the solution of the inverse prob¬ 
lem can be resolved from the measured data, and to 
which it has to be determined by the regularization, 
depends on the information content of the data [3]. If 
the information content is small, the effective number 
of degrees of freedom has to made small by strong 
regularization. If the information content is large, the 
amount of regularization can be diminished, allowing 
more discrimination in the solution. 

Adding more channels to an ECG or an MCG data set 
doesn’t increase the information content after a cer¬ 
tain point [4], as the new channels can then be de¬ 
scribed as a linear combination of the existing ones. 
Another possibility to increase the information con¬ 
tent of the data might be combining ECG and MCG 
data. The ECG and MCG are (in the UDL model) 
generated by the same source distribution, but as the 
sensitivity distribution over the heart differs for the 
ECG and MCG, one might expect that the informa¬ 
tion in the ECG and MCG is, to some extend, inde¬ 
pendent. 

In this study we investigate whether the results of the 
UDL-based inverse are improved by combining ECG 
and MCG data. The results of the UDL-based in¬ 
verse procedure obtained by using ECG and MCG 
data separately and combined were compared for a 
patient group for whom invasive data were available. 
Furthermore, in order to determine to what extend the 
MCG and ECG represent independent views of the 




Figure 1: Cross section of volume conductor model. 
S\: torso; S 2 , S 3 : lungs; £ 4 , S 3 : blood filled cavi¬ 
ties; Sh: ventricular myocardium. 


source distribution the transfer matrix of a combined 
ECG/MCG data set was compared to those of the sep¬ 
arate data sets. Finally, the information content in 
measured combined ECG/MCG data was compared 
to that in the separate measurement sets. 

2 Methods 

2.1 Patient data 

The study included 4 patients that had previously suf¬ 
fered myocardial infarction and that where scheduled 
for open chest surgery in order to treat ventricular ar¬ 
rhythmia. 

Prior to surgery, the ECG was recorded at 123 elec¬ 
trodes. The signals where sampled at 1000 Hz. Based 
on MRI images, an individual triangulated volume 
conductor model was constructed for each patient, in¬ 
cluding the torso (conductivity: 0.2 S/m), lungs (0.04 
S/m) and cavities within the heart (0.6 S/m). The 
ventricular surfaces, comprising both epi- and endo¬ 
cardium, were also included in the MRI based geom¬ 
etry model (see Fig. 1). 

During surgery, the potentials at the epicardium were 
recorded at 102 bipolar leads in an electrode sock 
system. The distance between the electrodes of the 
bipolar leads was 4 mm. Electrograms were mea¬ 
sured both for sinus beats and ectopic beats. In this 
study only the sinus beats were used. The electro¬ 
grams were sampled at 1000 Hz. From the electro¬ 
grams the activation times at the sites of the bipolar 
electrodes were determined. For monophasic electro¬ 
grams the activation time was defined as the time in¬ 
stant at which the absolute value of the potential was 
maximal, for biphasic electrograms it was defined as 




Figure 2: a) Early ventricular activation. Red: depo¬ 
larized myocardium. Blue: actual double layer. The 
arrows indicate the direction of the double layer, b) 
Equivalent double layer at ventricular surface, c) Ac¬ 
tual double layer at a later stage, d) Equivalent dou¬ 
ble layer at this stage. 


the time instant at which the absolute value of the 
derivative of the potential was maximal. Whether an 
electrogram was monophasic or biphasic was deter¬ 
mined visually. 

The locations of the epicardial electrodes within the 
geometry model were constructed by using an algo¬ 
rithm that simulates the positioning of the electrode 
sock around the heart on basis of a sketch made dur¬ 
ing surgery. The distribution of the activation times 
over the epicardium was estimated from the activation 
times at the epicardial electrodes by interpolation [5]. 

2.2 UDL source model 

During ventricular depolarization current is generated 
at the depolarization wavefront. Electrophysiology 
has demonstrated that the current source can be de¬ 
scribed as a uniform double layer with a strength of 
40 mV, located at this front [ 6 ]. 

Any uniform double layer with the same rim as the ac¬ 
tual wave front will generate the same external fields. 
Hence, the geometry of the wave front cannot be de¬ 
termined from measurements of these fields. In the 
UDL model an equivalent double layer is postulated 
at the ventricular surface Sh . This equivalent double 
layer comprises the region of the ventricular surface 
that has been reached by the depolarization wave front 




2.4 Combined inverse procedure 


(see Fig. 2). 

As can be seen in Fig. 2, during ventricular depolar¬ 
ization the equivalent source is completely described 
by activation function t(x) (which specifies the time 
at which activation reaches point x at Sh)- From that 
moment onward a 40 mV double layer is active at x. 

2.3 Electric and magnetic inverse procedure 

By using the boundary element method (BEM) the 
potential afj at observation point i at the body sur¬ 
face generated a uniform double layer with a strength 
of 40 mV at element j at the ventricular surface was 
computed (see figure 1 ), resulting in an electric trans¬ 
fer matrix A E . The BEM was also used to compute 
the magnetic transfer matrix A M {af- is the field gen¬ 
erated in gradiometer i by source element j ). For each 
patient the transfer matrices were computed using the 
patient’s individual volume conductor model. 

For the UDL source model the electric potential ip at 
electrode position i at time t is given by [7] 

N h 

¥>»(*) = ^2 a ij - T j) ) (!) 
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where Tj is the activation time at element j of Sh, and 
Nh is the number of elements of Sh- Similarly, the 
field at gradiometer i is given by 

N h 

B i {t) = Y J < j H{t-r j ). (2) 

3 = 1 

The electric inverse problem involves the estimation 
of the activation times r from the measured potentials 
cp. This problem is non-linear and ill-posed. It was 
solved by minimizing (by means of a quasi-Newton 
method) over all time samples simultaneously 

+Areg(-r), (3) 

where the regularization function reg(r) is defined as 

reg(r) = (^ (A(r(x ))) 2 dsfj 

(At denotes the surface Laplacian of r). The value 
of A was chosen such that reg(r) ~ 6 s/m. This 
corresponds to a physiologically realistic value for the 
velocity of the activation wave front of 0.17 m/s [7]. 
The initial estimate for r was obtained by applying 
the critical point theorem [ 8 ]. 

The equation for the magnetic inverse procedure is 
equal to equation 3, with ip replaced by B and A E by 
A M . 


When using the combined ECG and MCG data set 
in the inverse procedure, the data vectors (p (contain¬ 
ing the ECGs at all electrodes) and B (containing the 
MCGs at all gradiometers) have to be combined to 
one data vector C. Similarly, the transfer matrices 
have to be combined to one matrix: 




The scaling factor 7 is necessary because the ECG 
and MCG are of different dimension. 

The value of 7 follows from the formulation of 
the maximum likelihood estimation. Consider mea¬ 
surements yi with white noise of strength cr^, i = 
1, • • •, N. The maximum likelihood estimator of pa¬ 
rameter p in model yi = yi(p) is obtained by mini¬ 
mizing 

m - y%(p) 

If the noise level is the same for all measurements, the 
individual noise levels can be ignored. This is in fact 
the basis of equation 3, where the noise is assumed to 
be the same in all leads. 

From equation 4 can be seen that if the ECG and 
MCG channels contain white noise of strength cr E and 
a M respectively, the maximum likelihood solution is 
obtained if 7 = cr E /cr M . 

The equation for the combined inverse procedure is 
equal to equation 3, with cp replaced by C and A E by 
A c . 




2.5 Richness of the transfer matrices 

As the transfer matrix is not singular, the number of 
independent source patterns that can be distinguished 
at the heart is — in theory — equal to the number 
of channels in the measurement. However, many of 
these source patterns correspond to measurement sets 
that differ less than noise level. Hence only those k in¬ 
dependent patterns that correspond to measurements 
sets that differ more than the noise level can be dis¬ 
guised by the transfer matrix, k is called the effective 
rank of the transfer matrix, and depends on the noise 
level in the data, as well as the number and locations 
of recording sites. 

If, for two different recording set-ups, the effective 
rank of one of the corresponding transfer matrices is 
larger, this recording set-up can distinguish more in¬ 
dependent source patterns. In this paper, we will call 
this transfer matrix “richer” than the other one. 




Figure 3: Volume conductor model for patient 2. 
Some of the triangles at the front are not plotted to 
allow a view of the heart and the lungs. 


The singular value decomposition (S VD) of the trans¬ 
fer matrices provides a means to compare the rich¬ 
ness of transfer matrices. The SVD of a matrix A 
leads to the decomposition A = UYy 1 , where U and 
V are orthonormal matrices and E is a diagonal ma¬ 
trix. U and V are bases of the observation and source 
space, respectively. The diagonal elements cr* of E 
(also known as singular values) are a measure for the 
strength of the corresponding basis vectors. They are 
sorted such that 

The relative representation error r(k) is the relative 
difference between the actual transfer matrix, and the 
one in which all singular values beyond the effective 
rank (o^, i > k) are replaced by zero. If the represen¬ 
tation error r(k) of a matrix falls off fast as a function 
of k, the transfer matrix is predominantly represented 
by just a few basis vectors. Consequently, the source 
can generate only a few different observation patterns 
in this recording set-up. As a result, the inverse pro¬ 
cedure for such a transfer matrix is strongly ill-posed. 
In this study, we have used the speed with which r(k) 
falls off as a measure for the richness of the transfer 
matrix. 

When comparing transfer matrices, the one for which 
r(k) falls off slowest is the richer one, and con¬ 
sequently one can expect the corresponding inverse 
problem to be less ill-posed. 



Figure 4: Reconstructed positions of the epicardial 
electrodes for patient 2. Left: frontal view; right: pos¬ 
terior view. Only the electrodes with sufficient signal 
quality are shown. 


2.6 Information content 

Consider a set of m^ signals: m\(f)... m^(t). If 
this set contains k independent signals si(t)... Sk(t ), 
each of the mi can be described as a linear combina¬ 
tion of these independent signals plus noise: 

k 

m,i(t ) = ^2 hij Sj(t) + rii(t) (5) 

i=i 


The number of independent signals in a data set is 
called the information content of the data set [9]. 

The determination of the information content of a set 
of recorded signals involves the estimation of the min¬ 
imum number of independent signals needed to de¬ 
scribe this set. For this purpose, a singular value 
decomposition of the data matrix M (the rows of 
which are the signals) is performed. Let M(k) be 
the truncation of M to the first k singular values. If 
noise is white with level a the likelihood L(k) that 
M(k) = M is given by 


(' m ij - rhij(k )) 2 




where m is the number of time samples. 

The maximization of L(k) (or, equivalently the mini¬ 
mization of — log L(k)), does not lead to a correct es¬ 
timate of the information content, as obviously L(k) 
is maximal for k = £ (then Mk = Ml = M). 
To overcome this problem, the information criterion 
IC(k ) has been introduced [10], by adding a penalty 
term to the log-likelihood, correcting for the increased 
number of degrees of freedom for increasing values of 
k: 


IC(k,c n ) = -logL(fc) + c n u(k,£) , 

with c n a constant depending of the number of obser¬ 
vations n (= m • £), and v(k,l) the number of de¬ 
grees of freedom. It has been shown that v(k,l) = 
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Figure 5: Activation times for patient 2. Top row: anterior aspect, bottom row: posterior aspect. From left to 
right: measured activation times, activation times estimated from measured ECGs, from measured MCGs, and 
from the combined data set. Isochrones are drawn at 5 ms intervals. The circles indicate the positions of the 
epicardial electrodes. 


k{2£ — k + l)/2 [9] and that a robust estimator for 
the information content is obtained for the choice 
c n — \ logn [11]. This estimator of the number of 
independent signals in the data is called the Minimum 
Descriptor Length (MDL). In this study, the MDL 
was used to determine the information content in the 
ECG, MCG and combined data sets. 

3 Results 

3.1 Patient data 

Figure 3 shows the geometry model that was con¬ 
structed for patient 2 from the recorded MR images. 
Figure 4 shows the reconstructed positions of the epi¬ 
cardial electrodes for the same patient. Only those 
epicardial electrodes are shown for which the sig¬ 
nal quality was sufficient to allow the local activation 
times to be determined. Typically, about half of the 
102 epicardial electrodes had to be rejected because 
of poor signal quality. 

3.2 Estimated activation sequences 

In Fig. 5 the measured activation sequences for pa¬ 
tient 2 are depicted, as well as the estimated activa¬ 
tion sequences found from, respectively, the ECG, the 


MCG and the combined data set. The measured acti¬ 
vation times are determined by interpolation between 
the values at the epicardial electrodes, and are only 
plotted in regions close to these electrodes. 

Table 1 lists the RMS differences between the mea¬ 
sured and estimated activation times for all patients. 


Table 1: The RMS differences between the measured 
activation times and those estimated from the mea¬ 
sured ECGs, MCGs and from the combined data set. 


Patient 

Difference [ms] 
from ECG from MCG 

combined 

i 

27.7 

25.7 

24.9 

2 

14.0 

16.8 

12.9 

3 

11.3 

14.0 

10.0 

4 

27.8 

19.0 

24.4 

Average 

20.2 

18.9 

18.0 


3.3 Richness of transfer matrix 

In Fig. 6 the representation error r(k) is plotted on 
a log-scale as a function of the effective rank for the 
ECG, MCG and combined transfer matrices. From 




















Figure 6: Representation error of ECG (red, solid), 
MCG (green, dash-dotted) and combined (blue, 
dashed) transfer matrices as a function of the effec¬ 
tive rank. 


this figure it is clear that the representation error of the 
MCG transfer matrix error falls off much faster with 
the effective rank than the representation error for the 
ECG transfer matrix, indicating that the inverse pro¬ 
cedure based on MCG data is more ill-posed than the 
one based on ECG data. Also, the representation er¬ 
ror falls off slower for the combined transfer matrix 
than for the ECG transfer matrix, indicating that the 
combined transfer matrix is richer. 

3.4 Information content 

The number of independent signals in the measured 
data was estimated by using the MDL. The signal-to- 
noise ratio was estimated to be 1% for both the ECG 
and the MCG data. In table 2 the results are listed for 

Table 2: The number of independent signals in the 
measured data (as estimated by the MDL) for the 
ECG, MCG and combined data set. 


Patient 

(ecg) 

MDL 

(meg) (combined) 

i 

11 

8 

13 

2 

7 

5 

11 

3 

9 

8 

12 

4 

7 

10 

11 

Average 

8.5 

7.7 

11.8 


the ECG, MCG and combined data. 

4 Discussion 

As the positions of the epicardial electrodes were 
based on a sketch made during surgery, the locations 
are not very accurate. Also, as the electrograms were 
recorded in bipolar leads, the accuracy of the local 
activation times based on these electrograms is lim¬ 
ited. These factors reduce the value of the quantita¬ 
tive comparison between the measured and estimated 
activation times. Hence, the quantitative results listed 
in table 1 should be interpreted with some care. How¬ 
ever, these results do support the suggestion that for 
the combined data set the estimated activation se¬ 
quence corresponds better to the measured one that 
for the separate data sets. 

Qualitatively, there is a fairly good correspondence 
between the measured and estimated activation se¬ 
quences for patient 2 as displayed in Fig. 3. The mea¬ 
sured activation sequence shows two break-throughs 
at the anterior part of the epicardium. In the activation 
sequence based on the ECGs, also two break-throughs 
are present at the anterior epicardium, but at some¬ 
what different locations. This difference in locations 
might, for a large part, well the result of an error in the 
reconstructed position of the epicardial electrodes. In 
the ECG based estimated activation sequence a large 
late region is observed posteriorly, where the mea¬ 
sured activation is not so late. In this area the my¬ 
ocardium was damaged by infarction, and probably 
did not have normal activation. As the UDL model 
requires every part of the epicardial surface to be ac¬ 
tivated, the inverse procedure will find late activation 
times for regions where in reality the activation wave 
front has a reduced source strength. 

The activation pattern found from the MCG data 
set displayed in this figure is similar to the pattern 
found from the ECG. The main difference is that one 
break-though is located more posterior, and that there 
is a spot of late activation between the two break¬ 
throughs. 

In the activation sequence estimated from the com¬ 
bined data set there are again two anterior break¬ 
throughs. The late region at the anterior part of the 
outflow region of the left ventricle, which is visible 
in the measured activation sequence, is better repro¬ 
duced in the estimation based on the combined data 
set than in the ECG- or MCG-based estimation. 

The results for patient 2 as seen in Fig. 3 are typical 
in the sense that for all patients there is a fairly good 








Figure 7: Representation error of ECG (red, solid), 
original MCG (green, dash-dotted) and modified 
MCG (blue, dashed) transfer matrices as a function of 
the effective rank. The modified transfer matrix cor¬ 
responds to sensors that are closer to the surface. 


qualitative correspondence between the measured and 
estimated activation sequence, and that the qualitative 
correspondence is somewhat better for the estimation 
based on the combined data set. 

Fig. 6 shows that from an effective rank of about 10 
the representation error falls off slower for the com¬ 
bined ECG/MCG transfer matrix than for the ECG 
transfer matrix. The actual value of the effective rank 
depends on the signal-to-noise ratio in the data in a 
non-straightforward manner. However, as a first order 
estimate one can assume that the rank should at least 
be large enough to obtain a representation error that 
is smaller than or equal to the signal-to-noise ratio. 
For a signal-to-noise ratio of 1%, Fig. 6 yields effec¬ 
tive ranks of about 25 and 35 for the ECG and com¬ 
bined data set respectively. Hence, for realistic values 
of the signal-to-noise ratio the combined transfer ma¬ 
trix is richer (has a higher rank; can distinguish more 
independent source patterns) than the ECG transfer 
matrix. 

Fig. 6 also shows that the representation error for the 
MCG transfer matrix falls off much faster than that of 
the ECG transfer matrix. This is not the result of some 
intrinsic inferiority of the MCG, but of the fact that 
the MCG sensors are further away from the source 
than the ECG sensors. Any source distribution can 
be described as the sum of a dipole, quadrupole, etc. 
The field generated by the dipole components falls off 
with the square distance, the quadrupole field with the 
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Figure 8: Representation error of ECG (red, solid), 
original MCG (green, dash-dotted) and modified 
MCG (blue, dashed) transfer matrices as a function 
of the effective rank. The set-up corresponding to the 
modified transfer matrix contains coils oriented in all 
three directions, positioned at the same locations as 
in the original set-up. 


cube of the distance, etc. So the larger the distance 
from the sensors to the source, the more the measured 
field will be dominated by the (3 components of the) 
dipole component, and the less rich the corresponding 
transfer matrix will be. 

To verify this we have constructed an MCG transfer 
matrix for a sensor set-up where the coils are 3.5 cm 
closer to the thorax than in the actual MCG set-up. 
The results are shown in Fig 7. The representation 
error for the MCG sensors that are closer to the tho¬ 
rax falls off much slower that for the original MCG 
sensors, but still faster than for the ECG set-up. This 
results from the fact that for the new sensor locations 
the sensors are still, in average, further away from the 
source than the ECG sensors, and also from the fact 
that the MCG sensor locations have a limited spatial 
spread. 

We have also investigated the influence of the coil ori¬ 
entations on the richness of the transfer matrix. In the 
original set-up, most coils were oriented more or less 
parallel to the body surface. We have constructed an¬ 
other observation set-up where at each location coils 
were added that were oriented in the two planes per¬ 
pendicular to the original coil. The results are plot¬ 
ted in Fig. 8. This figure shows that, by adding these 
perpendicular coils, the richness of the corresponding 
transfer matrix increases. 




The results in table 2 show that the information con¬ 
tent of the data increases when the ECG and MCG 
data are combined. To test whether this is an expres¬ 
sion of some kind of independence of the ECG and 
MCG data, or just an effect of the larger number of 
channels in the combined data set, the information 
content for patient 2 was also computed for a reduced 
ECG data set. The reduced data set (containing 60 
leads), was constructed by removing all but one odd- 
numbered leads from the complete data set (contain¬ 
ing 118 leads). The results are listed in table 3. 

Table 3: The influence of the number of leads on the 
number of independent signals for patient 2. The top 
3 rows repeat the results listed in table 3 for this pa¬ 
tient, and specifies the number of channels. The bot¬ 
tom 2 rows show the results when the ECG data set is 
replaced by a reduced ECG data set. 


data set 

# channels 

MDL 

ecg 

118 


7 

meg 

58 


5 

eeg+meg 

176 

(118+58) 

11 

ecg' 

60 


7 

eeg'+meg 

118 

(60+58) 

10 


This clearly shows that from a certain number of 
channels onward adding more leads to an ECG data 
set does not yield more information about the source. 
Also, it shows that adding MCG channels to an ECG 
data set does increase the information content of the 
data. 

Overall, we conclude that the transfer matrix for a 
combined ECG/MCG measurement set-up is larger 
than those for the ECG and MCG separately. Con¬ 
sequently, we have found in measured data that the 
information content in a combined ECG/MCG data 
set is larger than in each separately. Hence, the com¬ 
bination of ECG and MCG data in inverse procedures 
will improve the performance. The patient data pre¬ 
sented in this paper do not allow a precise quantitative 
assessment, but both the qualitative and quantitative 
evaluation do support this conclusion. 
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